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Abstract 

In this paper structure-preserving time-integrators for rigid body-type me- 
chanical systems are derived from a discrete Hamilton-Pontryagin variational 
principle. From this principle one can derive a novel class of variational parti- 
tioned Runge-Kutta methods on Lie groups. Included among these integrators 
are generalizations of symplectic Euler and Stormer-Verlet integrators from flat 
spaces to Lie groups. Because of their variational design, these integrators pre- 
serve a discrete momentum map (in the presence of symmetry) and a symplectic 
form. 

In a companion paper, we perform a numerical analysis of these methods 
and report on numerical experiments on the rigid body and chaotic dynamics of 
an underwater vehicle. The numerics reveal that these variational integrators 
possess structure-preserving properties that methods designed to preserve mo- 
mentum (using the coadjoint action of the Lie group) and energy (for example, 
by projection) lack. 
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1 Introduction 

Overview. This paper is concerned with efficient, structure-preserving time inte- 
grators for mechanical systems whose configuration space is a Lie group based on the 
Hamilton-Pontryagin (HP) variational principle Livens [1919]; Lall & West [2006]; 
Kharevych et al. [2006]; Yoshimura and Marsden [2006a, b]. This HP Principle 
has many attractive theoretical properties; for instance, how it handles degenerate 
Lagrangian systems. The present paper paper shows that the HP viewpoint also 
provides a practical way to design discrete Lagrangians, which are the cornerstone of 
variational integration theory. This overview explains the central idea of this paper 
in the context of vector spaces and shows how this idea extends to Lie groups. 

The HP principle states that a mechanical system traverses a path that extrem- 
izes the following HP action integral: 

&HP= / L{q,v)dt+ / {p,q-v)dt. (1.1) 

J a J a 

Lagrangian kinematic constraint 

The integrand of the HP action integral consists of two terms: the Lagrangian and a 
kinematic constraint paired with a Lagrange multiplier (the momentum). The kine- 
matic constraint relates the mechanical system's velocity to a curve on the tangent 
bundle. In this principle, the curves g(t), v{t\ p{t) are all varied independently. If 
p is varied first, it collapses to the usual Hamilton principle. If, on the other hand, 
v(t) is varied first it defines the (negative of the) Hamiltonian as the extrema of the 
terms involving v and then the principle reduces to Hamilton's phase space principle. 
This HP form of the action integral makes it amenable to discretization. 

In particular, one can implement an s-stage Runge-Kutta (RK) discretization of 
the kinematic constraint and enforce this discretization as a constraint in a discrete 
action sum. The motivation is that the theory, order conditions, and implementation 
of such methods, are mature. For this purpose let [a, b] and N be given, and define 
the fixed step size h = (6 — o) /N and tk = hk + a, k = 0, N . Let s be the number 
of stages in the RK method. In analogy with the continuous system, the discrete 
HP action sum takes the following form: 

N~l s 



k=0 i=l 



discrete Lagrangian 
N-1 s 



fc=0 i=l \ j=l I \ j=l I 



discrete kinematic constraint 

(1.2) 



It consists of two parts: a weighted sum of the Lagrangian using the weights from 
the Butcher tableau of the RK scheme, and pairings between discrete internal and 
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external stage Lagrange multipliers and the discretized kinematic constraint. This 
strategy yields a Lagrangian analog of a well-known class of symplectic partitioned 
Runge-Kutta methods including the Lobatto IIIA-IIIB pair which generalize to 
higher-order accuracy Suris [1990]; Marsden & West [2001]; Hairer, Lubich, and 
Wanner [2006]. 

In the Lie group context, one can generalize this strategy using either constrained 
or generalized coordinates. To use constrained coordinates one treats the system as 
a holonomically constrained mechanical system. In this approach one assumes that 
G can be written as the level set of some function g : M" M.^ , embeds G in a larger 
linear space, and uses Lagrange multipliers to enforce the constraint. This approach 
is discussed in Bou-Rabee & Owhadi [2007b]. The corresponding constrained action 
takes the following form: 



fc=0 i=l 



h 



(1.3) 



In the present paper a second approach based on generalized coordinates is 
presented. First the paper introduces the following left-trivialized action: 

5hP= ['i{g,Odt + f\ij,g-^g-0dt. (1.4) 

J a J a 

left-trivialized Lagrangian reconstruction equation 

Then an equivalence is established between critical points of Sjfp and &hp- If 
the Lagrangian is left-invariant, it is shown that this principle unifies the system's 
Lie-Poisson and Euler-Poincare descriptions Marsden & Scheurle [1993]; Cendra, 
Marsden, Pekarsky, and Ratiu [2003]. Since the reconstruction equation is a differ- 
ential equation on a Lie group, one cannot directly discretize it by an RK method. 
However, one can discretize it using an s-stage Runge-Kutta-Munthe-Kaas (RKMK) 
method Munthe-Kaas [1995]; Munthe-Kaas & Zanna [1997]; Munthe-Kaas [1998]; 
Munthe-Kaas & Owren [1999]. The integral of the left-trivialized Lagrangian is 
approximated using a weighted sum given by the 6-vector in the Butcher tableau 
of the RKMK scheme. This approach is shown to yield a novel class of variational 
partitioned Runge-Kutta (VPRK) methods on Lie groups; including generalizations 
of symplectic Euler and Stormer-Verlet methods on flat spaces. 



2 Background and Setting 

In the next paragraphs we will give some background material for the reader's 
convenience as well as to put the paper into context. 

Variational Integrators. Variational integration theory derives integrators for 
mechanical systems from discrete variational principles. The theory includes discrete 
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analogs of the Lagrangian, Noether's theorem, the Euler-Lagrange equations, and 
the Legendre transform. Variational integrators can readily incorporate holonomic 
constraints (via Lagrange multipliers or the discrete null-space method; Leyen- 
decker, Marsden, and Ortiz [2007]) and non-conservative effects (via their virtual 
work) Marsden & West [2001], as well as discrete optimal control (see Leyendecker, 
Ober-Blobaum, Marsden, and Ortiz [2007] and references therein). Altogether, this 
description of mechanics stands as a self-contained theory of mechanics akin to 
Hamiltonian, Lagrangian or Newtonian mechanics. 

One of the distinguishing features of variational integrators is their ability to 
compute statistical properties of mechanical systems, such as in computing Poincare 
sections, the instantaneous temperature of a system, etc. For example, as a conse- 
quence of their variational design, variational integrators are symplectic. A single- 
step integrator applied to a mechanical system is called symplectic if the discrete 
flow map it defines exactly preserves the canonical symplectic 2-form and is other- 
wise called standard. Using backward error analysis one can show that symplectic 
integrators applied to Hamiltonian systems nearly preserve the energy of the contin- 
uous mechanical system for exponentially long periods of time and that the modified 
equations are also Hamiltonian Hairer, Lubich, and Wanner [2006]. Standard inte- 
grators often introduce spurious dynamics in long-time simulations, e.g., artificially 
corrupt chaotic invariant sets is well-illustrated in a computation from Bou-Rabee 
[2007], namely of a Poincare section of an underwater vehicle in Fig. 2.1 using a 
fourth-order accurate Runge-Kutta (RK4) method and a variational Euler (VE) 
method designed for rigid-body type systems. 




Figure 2.1: Underwater Vehicle Dynamics. This figure shows a computation of 
Poincare sections using a second-order accurate variational Euler integrator (VE) as com- 
pared to fourth-order accurate Runge-Kutta (RK4). Both methods agree with the bench- 
mark at the finer stepsize h = 0.025. However, at the coarser stepsize h = 0.05, RK4 
corrupts chaotic invariant sets while the lower-order accurate VE method preserves the 
structure of the benchmark. 



In addition to correctly computing chaotic invariant sets and long-time excellent 
energy behavior, evidence is mounting that variational integrators correctly compute 
other statistical quantities in long-time simulations. For example, in a simulation of 
a coupled spring-mass lattice. Lew, Marsden, Ortiz, and West [2004] found that vari- 
ational integrators correctly compute the time-averaged instantaneous temperature 
(mean kinetic energy over all particles) over long-time intervals, whereas standard 
methods (even a higher-order accurate one) exhibit a artificial drift in this statistical 
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quantity. These structure- preserving properties of variational integrators motivated 
their extension to stochastic Hamiltonian systems. 

Structure-Preserving Lie Group Integrators. For a mechanical system on a 
Lie group that possesses the symmetry of that Lie group, in addition to the symplec- 
tic structure, the resulting flow preserves a momentum map associated with the Lie 
group symmetry. In this context there are several different strategies available to 
derive structure-preserving Lie group integrators; some of these are discussed here. 

One strategy involves the so-called Lie-Newmark method due to Simo & Vu- 
Quoc [1988] and Simo & Wong [1991]. These methods were motivated by the need 
to develop conserving algorithms that efficiently simulate the structural dynamics of 
rods and shells. For example, the configuration space of a discrete, three-dimensional 
finite-strain rod model, would involve N copies of x S0(3) where is the number 
of points in the discretization of the line of centroids of the rod. For each point on 
the line of centroids, the orientation of the rod at that point is specified by an 
element of S0(3). In such models the mathematical description of the rotational 
degrees of freedom at these points is equivalent to the EP description of a free rigid 
body with added nonconservative effects due to the elastic coupling between points. 

It was not apparent that the proposed Lie-Newmark methods had the necessary 
structure-preserving properties. In fact, Simo & Wong proposed another set of algo- 
rithms which preserve momentum by using the coadjoint action on S0(3) to advance 
the flow. Such integrators will be referred to as coadjoint-preserving methods. Later, 
Austin et al. [1993] showed that the midpoint rule member of the Lie-Newmark 
family with a Cayley reconstruction procedure was, in fact, a coadjoint-preserving 
method for S0(3). They also numerically demonstrated the method's good per- 
formance crediting it to third-order accuracy in the discrete approximation to the 
Lie-Poisson structure. In related work, McLachlan & Scovel [1995] construct re- 
duced, coadjoint-orbit preserving integrators by reducing G-equivariant integrators 
on T*G obtained by embedding G in a linear space using holonomic constraints. 

Coadjoint and energy preserving methods of the Simo & Wong type that fur- 
ther preserve the symplectic structure were developed for S0(3) by Lewis & Simo 
[1994, 1996]. This was done by defining a one-parameter family of coadjoint and 
energy-preserving algorithms of the Simo & Wong type in which the free parame- 
ter is a functional. The function was specified so that the resulting map defined a 
transformation which preserves the continuous symplectic form. 

Endowing coadjoint methods with energy-preserving properties was also the sub- 
ject of Eng0 Sz Faltinsen [2001]. Specifically, they introduced integrators of the 
Runge-Kutta Munthe-Kaas type that preserved coadjoint orbits and energy using 
the coadjoint action on S0(3) and a numerical estimate of the gradient of the Hamil- 
tonian. 

Variational integration techniques have been used to derive structure-preserving 
integrators on Lie groups; see Moser & Veselov [1991]; Wendlandt and Marsden 
[1997]; Marsden, Pekarsky, and Shkoller [1998]; Bobenko and Suris [1999a]; Bobenko 
& Suris [1999b]. Moser and Veselov derived a variational integrator for the free 
rigid body by embedding S0(3) in the linear space of 3 x 3 matrices, M^, and using 
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Lagrange multipliers to constrain the matrices to S0(3). This procedure was subse- 
quently generalized to Lagrangian systems on more general configuration manifolds 
by the introduction of a discrete Hamilton's principle on the larger linear space with 
holonomic constraints to constrain to the configuration manifold in Wendlandt and 
Marsden [1997]. They also considered the specific example of deriving a variational 
integrator for the free rigid body on the Lie group by embedding into and 
using a holonomic constraint. The constraint ensured that the configuration update 
remained on the space of unit quaternions (a Lie group) and was enforced using a 
Lagrange multiplier. 

Another approach is to use reduction to derive variational integrators on re- 
duced spaces. Marsden, Pekarsky, and ShkoUcr [1998] developed a discrete analog 
of EP reduction theory from which one could design reduced numerical algorithms. 
They did this by constructing a discrete Lagrangian on G x G that inherited the 
G-symmetry of the continuous Lagrangian, and restricting it to the reduced space 
(G X G)/G ~ G. Using this discrete reduced Lagrangian and a discrete EP (DEP) 
principle, they derived DEP algorithms on the discrete reduced space. They also 
considered using generalized coordinates to parametrize this discrete reduced space, 
specifically the exponential map from the Lie algebra to the Lie group. These tech- 
niques were applied to bodies with attitude-dependent potentials, discrete optimal 
control of rigid bodies, and to higher-order accuracy in Leok, McClamroch , and 
Lee [2005] and Lee, Leok, and McClamroch [2007]. 

Bobenko and Suris [1999a] considered a more general case where the symmetry 
group is a subgroup of the Lie group G in the context of semidirect Euler-Poincare 
theory (see Holm et al. [1998]). They did this by writing down the discrete Euler 
Lagrange equations for this system and left-trivializing them. For the case when 
the symmetry group is G itself, one recovers the DEP algorithm as pointed out in 
Marsden, Pekarsky, and Shkoller [1998]. In addition, Bobenko & Suris [1999b] used 
this theory to determine and analyze an elegant, integrable discretization of the 
Lagrange top. 

The perspective in this paper on Lie group variational integrators is different. 
Recognizing that Euler's equations for a rigid body are in fact decoupled from the 
dynamics on the Lie group, and more generally, that the EP equation is decoupled 
from the dynamics on the Lie group, the paper aims to develop discrete variational 
schemes that analogously consist of a reconstruction rule and discrete EP equations 
that can be solved independently of the reconstruction equation and on a lower di- 
mensional linear space. As mentioned in the overview the central idea is to discretize 
the reduced HP principle. 

Organization of the Paper. In §3 continuous HP mechanics and its reduction is 
presented. In particular, it is shown that the reduced and unreduced HP variational 
principle are equivalent to Hamilton's and the EP variational principles. Moreover, 
properties of the HP flow map are verified mainly to guide the discrete theory. In 
§4 the reduced discrete analog of the HP theory is developed. Properties of the 
discrete flow map are verified including discrete momentum map and symplectic 
form preservation. The theory is illustrated on several specific examples. In §5 the 
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structure-preserving Lie group integrators relevant to this paper are presented. In 
§6 the free rigid body and underwater vehicle examples are presented, the structure- 
preserving methods from §5 specialized to these examples, and results of numerical 
experiments are presented. 

Part II of this Paper. The second installment of this paper will be devoted to 
the numerical analysis of HP methods along with numerical experiments on a class 
of nonreversible mechanical systems on Lie groups as well as the chaotic dynamics 
of an underwater vehicle. A specific outline of that paper is given in the conclusion 
section of the present paper. 



This section develops basic mechanics on Lie groups from the Hamilton- Pontryagin 
perspective. 

The HP Principle. Consider a mechanical system whose configuration space is 
a Lie group G. Let its tangent and cotangent bundles be denoted TG and T*G 
respectively, and its Lie algebra and dual be given by q and q* respectively. In this 
paragraph the left-trivialization of the HP principle for a Lagrangian L : TG M 
will be derived. 

The HP principle unifies the Hamiltonian and Lagrangian descriptions of a me- 
chanical system, as shown in Yoshimura and Marsden [2006a, b]. It states the fol- 
lowing critical point condition on TG © T*G, 



where {g(t),v{t),p{t)) £ TG ® T*G are varied arbitrarily and independently with 
endpoint conditions g{a) and g{b) fixed. This builds in the Legendre transformation 
as well as the Euler-Lagrange equations into one principle. 

Definition 3.1. Following standard conventions, the left action of G on TG orT*G 
is denoted by simple concatentation. The left-trivialized Lagrangian £ : G x g ^ M 
is defined as: 



The HP principle for mechanical systems on Lie groups is equivalent to the left 
trivialized HP principle: 
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l{g,i)=L{g,gO. 




where there are no constraints on the variations; that is, the curves ^(t) G g, ^(t) G 
g* and g{t) G G can be varied arbitrarily. To see this, we proceed as follows. 
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Let S{g,v,p) denote the HP action functional or integral, 
Sig,v,p)= f [L{g,v) + {p,g-v)]dt. 

J a 

Fixing the interval [a, b], we regard S as a map on path space: S : C{TG®T*G) — > M, 
where 

C{TG®T*G) = {{g,v,p) : [a,h]^TG®T*G \ {g,v,p) £ C°°{[a,b],TG ®T*G)}. 
Then a simple calculation shows that, 

rb 

S{g,v,p)= / [L{g,gO + {p,gg-\g-v))] dt 

J a 
J a 
J a 

= si9,^,lJ^) 



dt 
dt 



where s is the reduced HP action functional, = g ""^u G fl, and n = g G fl*. 
Prom this equality one can derive the following key theorem. 

Theorem 3.2. Consider a Lagrangian system on a Lie group G with Lagrangian 
L : TG M. Let £ : G x q ^ W be its left-trivialization. Then the following are 
equivalent 

1. Hamilton's principle for L on G 

5 [ Lig,g)dt = 

J a 

holds, for arbitrary variations g{t) with endpoint conditions g{a) and g{b) fixed; 

2. the following variational principle holds on Q, 

S f l{g,Odt = Q 

J a 

using variations of the form 

= ?7 + ad^ rj 

where rj{a) = r]{b) = and ^ = g~^g; i.e., $, = TLg-ig; 

3. the HP principle 

S [ [Lig,v) + {p,g-v)] dt = 

J a 

holds, where {g{t),v(t),p{t)) G TG ®T*G, can be varied arbitrarily and inde- 
pendently with endpoint conditions g{a) and g{b) fixed; 
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4- the left-trivialized HP principle 

rb 



J a 



holds, where {g{t),^{t), fi{t)) G G x g x g* can be varied arbitrarily and inde- 
pendently with endpoint conditions g{a) and g{h) fixed. 



Remark. If the Lagrangian is left-invariant, i.e., L(g,v) = L(hg, hv) for all h £ G, 
then the left-trivialized Lagrangian simplifies. In particular, taking h = g~^, = 
L{9i90 — -^(C)Oi where e is the identity element of the group. In this case the 
left-trivialized HP principle unifies the Euler-Poincare and Lie-Poisson descriptions 
on Q and g* respectively, consistent with the results of Marsden &; Scheurle [1993] 
and Cendra, Marsden, Pekarsky, and Ratiu [2003]. 

The HP Flow. From the left-trivialized HP principle, the variations of s with 
respect to and fi give 



varying /i gives ^ = 9 ^9, (reconstruction equation) 
d£ 

varying ^ gives = -Q^i9,0i (Legendre transform). 
Also, setting the variation of s with respect to g equal to zero gives 



d£ 
dg' 



,59) + {l^,5{9 g) 



dt 



9-Qg^9 ^59) + {fJ;-9 ^Sgg ^g + g ^6g) 



dt = 



(3.1) 
(3.2) 



(3.3) 



Observe that 



[{fi,6{g-^g))] dt= [{fi,-g-^6gg~^g + g-^5g)] dt 



Let r] = g ^5g. Using the product rule and (3.1), we see that 



dt 



-iv + 9 ^'rJd, which implies g ^'^5g = '^V + iV- 
at dt dt 



Substituting this relation into (3.3) gives 



dl \ / d 



dt = 0. 



Integration by parts and using the boundary conditions on g yields 



d 
dt 



+ ad^ H + g—,ri 



dt = 0. 
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Since the variations are arbitrary, one arrives at 



dt 



fj, = ad^ fi + g 



di 
dg' 



(3.4) 



In sum, the left-triviahzed HP equations are given by: 




ad| /X + 5 



di 



(3.5) 



Assuming that the Legendre transform is invertible, (3.5) describes an IVP on the 
left-triviahzed space G x g x g*. 

Definition 3.3. Letlnp denote the admissible space and defined as, 



Let Ihp denote its left-trivialization and defined as the subset of G x g x q* that 
satisfies (3.2), i.e., 



The natural projection is denoted by tthp '■ TG ®T*G ^ T*G and defined as, 

T^Hp{g,v,p) := {g,p), '^H^p{g,p) = [g,v,p), {g,v) =¥L'^{g,p) 
where ¥L is the Legendre transform. 

Given a time-interval [a, b] and an initial ((^(a), ^(a), /x(a)) G Ihp, one can solve for 
{g{b),(,{b), fi{b)) £ Ihp by eliminating ^ using the left-trivialized Legendre transform 
(3.2) and solving the ODEs (3.1) and (3.4) for g and fj,. Let this map on Ihp be 
called the left-trivialized HP flow map, Fhp : Ihp Ihp- 

The flow map Fhp is equivalent to the HP flow on Ih p through left trivialization 
which defines a diffeomorphism between TG ®T*G and G x g x g*, and hence, 
between Ihp and Ihp- Through tthp the HP flow is identical to the Hamiltonian 
flow for the Hamiltonian of this mechanical system on T*G obtained via the Legendre 
transformation. Although it hp is not a diffeomorphism from TG ®T*G to T*G, it 
is a diffeomorphism when its domain is restricted to Ihp- Thus, the left-trivialized 
HP, HP and Hamiltonian flows of this mechanical system are all equivalent. This 
observation makes the subsequent proof of symplecticity seem superfluous, since this 
structure obviously follows from the standard theory of Hamiltonian systems with 
symmetry. However, this verification is still important since it serves as a model for 
the less obvious discrete theory. 

It will be helpful to define 'KXhp = t^hp\ihp- The manifold TG ®T*G is a 
presymplectic manifold with the HP presymplectic form, ^hp = t^^p^-, and the 




(3.6) 




(3.7) 
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manifold Ihp is a symplectic manifold with the HP symplectic form, ^i^p = vtj^^O. 
Similarly, the manifold Gxqxq* is a presymplectic manifold with the presymplectic 
form ujhp that is obtained by puUing-back the HP presymplectic form by the left 
trivialization of TG T*G, <p : G x q x q* ^ TG (B T*G, i.e., lohp = 4>*^hp- 
However, if the left-trivialization is restricted to Ihp, 4'ihv ~ 4'\ihp^ then 2hp is a 
symplectic manifold with the symplectic form given by ujx^j, = (l^x^ ^Thp- 

Symplecticity. The symplectic structure of left-trivialized HP flows is obvious 
from the standard theory of Hamiltonian systems with symmetry, but reviewing the 
proof will help since it parallels the discrete case. 

Consider the restriction of the left-trivialized HP action integral to solutions of 
(3.5): s. Since the space of solutions of (3.5) can be identified with X^p, s : X^p M. 
The differential of s can be written as. 



ds • {5g{a),dC{a),6fj,{a)) 



dt 



+ 



dt+ {n,g'^6g) 



a 



= {fi,g~'5g)\l = mpr9j,^-6ij • {6g{a),SC{a),6fi{a)) 

where we have introduced the left-trivialized HP one-form, Ojf^^ = ©x^p- Since 
d^s = 0, observe that 

d2s = (F,p)*o;x,,-u;i,, = 0. 
And hence, as a map on X^p, F^p is symplectic. 

Theorem 3.4. Left-trivialized HP flows preserve the symplectic two-form coj^p- 

4 Lie Group VPRK Integrators 

The purpose of this section is to use the general HP methodology to derive a variety 
of integrators of variational partitioned Runge-Kutta (VPRK) type for Lie groups. 
After introducing the map r which is typified by the exponential map, and its 
properties, we use an s-stage Runge-Kutta-Munthe-Kaas (RKMK) approximation 
to the reconstruction equation, which leads naturally to the introduction of VPRK 
Integrators on Lie groups. This includes the Stormer-Verlet method for Lie groups, 
variational Euler methods on Lie groups, and Euler-Poincare integrators. 

Canonical Coordinates of the First Kind. To setup the discrete HP principle, 
we introduce a map r : g — > G. Let e G G be the identity element of the group. 
The map r is assumed to be a local diffeomorphism mapping a neighborhood of zero 
on g to one of e on G with t(0) = e, assumed to be analytic in this neighborhood, 
and assumed to satisfy r(^) • t(— ^) = e. Thereby r provides a local chart on the 
Lie group. By left translation this map can be used to construct an atlas on G. 
An example of a r is the exponential map on G, but there are other interesting 
examples as well, as we shall see shortly. 
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Definition 4.1. The local coordinates associated with the map t are called canon- 
ical coordinates of the first kind or just canonical coordinates. 

For an exposition of canonical coordinates of the first and second kind, and their 
apphcations the reader is referred to Iserles, Munthe-Kaas, N0rsett, and Zanna 
[2000] . In what follows we will prove some properties of these coordinates that will 
be needed shortly. 

Derivative of r and its inverse. To derive the integrator that comes from a 
discrete left-trivialized HP principle, we will need to differentiate t~^. The right 
trivialized tangent of r and its inverse will play an important role in writing this 
derivative in an efficient way. The following is taken from Definition 2.19 in Iserles, 
Munthe-Kaas, N0rsett, and Zanna [2000]. 

Definition 4.2. Given a local diffeomorphism t : q ^ G, we define its right 
trivialized tangent to be the function dr : g x g — > g which satisifies, 

DT{0-6 = TR,^^)dT^{6). 

The function dr is linear in its second argument. 

Figm'e 4.1 illustrates the geometry behind this definition. 




Figure 4.1: Derivative of r. Definition (4.2) splits tlie difTerential of r into a map on tlie Lie 
algebra (the right trivialized tangent of r) and right multiplication to the tangent space at r(^). 



From this definition the following lemma is deduced. 
Lemma 4.3. The following identity holds, 



dT^{6) = AdT-(g) (ir_g(5). 
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Proof. Differentiation of t(i^) • t(— ^) = e gives 

D r(-e) • 5 = -TL,^_^)TR,^^^) (D r(0 • S) . 
While tlie chain rule yields 

Dr(-e)-5 = -ri?,(_^)dT_5(5). 
Combining these two identities and using the definition above, 

Simplifying this expression gives, 

which proves the identity. ■ 
We will also need a simple expression for the differential of r^^. 

Definition 4.4. The inverse right trivialized tangent ofr is the function dr^^ : 
^ which satisifies for g = r(^), 

BT-^{g) • 6 = dT^^{TR^(^_i:)6), dT^^{dTi.{6)) = 6. 

The function dr^^ is always linear in its second argument. 
Figure 4.2 illustrates the geometry behind this definition. 




Figure 4.2: Derivative of r ^. Definition 4.4 splits tlie differential ofr ^ into riglit multiplica- 
tion to the Lie algebra and a map on the Lie algebra (the right trivialized tangent of r^^). 

The following lemma follows from this definition and Lemma 4.3 above. 
Lemma 4.5. The following identity holds, 

dT^\6) = drr| (Ad^(_g) 6). 



4 Lie Group VPRK Integrators 



14 



Proof. This follows directly from Lemma 4.3. Let 6 dT^^{6) in that identity to 
obtain 

And now solve this equation for dT^^{6), 

dT~\6) = drZl (Ad^(_^) 5) . 



RKMK Discretization of Reconstruction Equation. Let [a, b] and N be 

given, let h = {b — a)/N be a fixed integration time step and = hk. A good 
candidate for discretizing the reconstruction equation is given by a generalization 
of s-stage Runge-Kutta methods to differential equations on Lie groups, namely 
Runge-Kutta-Munthe-Kaas (RKMK) methods introduced in the following series of 
papers: Munthe-Kaas [1995]; Munthe-Kaas & Zanna [1997]; Munthe-Kaas [1998]; 
Munthe-Kaas & Owren [1999]. The idea behind those papers is to use canonical 
coordinates on the Lie group to transform the differential equation on TG, e.g., 
given by, 

g = gf{t,g), g{0) = go, g{t) e G, f{t,g{t))e9. (4.1) 

to a differential equation on g. Specifically, substitute the following parametrization 
9{t) = 5o'^(0(i)) into (4.1) to obtain, 

9 = TLg,;TR^(Q)dTeQ = TLg^,TL^^0)f{t, g). 

Using Lemma 4.3 this equation can be rewritten as, 

TL^(^_Q)TR^(^0)dTe& = Ad^(_e) drsQ = dr^eQ = f{t,g). 

Solving for Q gives 

e = dT:y{t,g), e(0)=0, e(t)Gg. (4.2) 

As described in the following definition, the RKMK method is obtained by applying 
an s-stage RK method to (4.2). 

Definition 4.6. Consider the first-order differential equation g = f{t, g) for the 
curve {g{t), f{t, g{t))) £ TG. Given coefficients bi,aij £ M (i,j = 1, ■ ■ ■ , s) and set 
Ci = "^j^iO-ij- An s-stage Runge-Kutta-Munthe-Kaas (RKMK) approximation 
is given by 

G{ = gkT{h@l), (4.3) 

s 

©1 = ^ E ^i^dT2l^J{tk + Cjh, Gi), i = 1, • • • , s, (4.4) 
9k+i = 5fcT f /i ^ bjdT^l^Jitk + cjh, G{) \ . (4.5) 
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If aij = for i < j the RKMK method is called explicit, and implicit otherwise. The 
vectors and G\ are called external and internal stage configurations, respectively. 



It follows that for given r an s-stage RKMK method is determined by its a-matrix 
and 6-vector which are typically displayed using the so-called Butcher tableau: 



Cl 


ail ■ 


■ ais 


Cs 


asi ■ 


O-ss 




hi ■ 


■ hs 



Suppose that ^(t), t S [a, 6], is given. From this definition it is clear that an s-stage 
RKMK method applied to g = g^ can be written as: 



(4.6) 



where = ^(t^ + Cih). In practice one often truncates the series expansion of 
dT~^^j . The following theorem guides how to do this without degrading the order 

of accuracy Hairer, Lubich, and Wanner [2006]. 

Theorem 4.7. Given a qth order approximant to the exact exponential: t : q ^ G. 
If the underlying RK method is of order p and the truncation index ofdrZQ satisfies 
q > p — 2 then the RKMK method is of order p. 

VPRK Integrators on Lie Groups. The discrete HP principle states that the 
discrete path the discrete system takes is one that extremizes a reduced action sum 
that will be introduced shortly. To discretize the action integral, (4.6) is treated 
as a constraint in the discrete HP action, and the integral of the left-trivialized 
Lagrangian is approximated by the following quadrature: 



rtk+h ^ 

/ £{g,Odt^Y.^b^^('Gl^^k)- (4.7) 
Jtk i^i 



The truncation index of dr_Q in (4.6) is chosen to he q = 0. By theorem 4.7 one 
can obtain second-order accurate methods from this principle. 

Definition 4.8. Given an s-stage RKMK method with bj ^ for j = 1, s, define 
the discrete VPRK path space, 

Cd{gi,g2) = {ig,fiAQ\^\l^'}t=i)d ■■ {tk}to - (G x g*) x (g x g x g*y \ 

g{to)=gi, g{tN) = g2}- 
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and the action sum Sd : Cd{gi-,g2) 

N-l s 



as 



k=0 1=1 



\ 3=1 I 



{^^k+i,r ^{g^^gk+i)/h-^bjE- 



(4.8) 



Observe that is an approximation of the reduced HP action integral by numer- 
ical quadrature. The definition of r as a map from g to G ensures that the pairings 
in the above sum arae well defined. The discrete left-trivialized HP principle states 
that, 

Ssd = 

for arbitrary and independent variations of the external stage vectors {gk, fJ-k) £ 
G X Q* and the internal stage vectors (6^, E\, G g x g x g* for i = 1, • • • , s and 
k = 0, - ■ ■ ,N subject to fixed endpoint conditions on {^fcj^o- 

Theorem 4.9. Let i : Gxg ^ M 6e a smooth, left-trivialized Lagrangian. A discrete 
curve Cd € Cd{gi,g2) satisfies the following VPRK scheme: 



r-\g-'Gi)/h = J2a,jEi = @i, 



s 



(4.9) 
(4.10) 



r ^(5fc ^gk+i)/h = bjEl = ^k+i, 

+ (^.('^^eV* - "Ifidr^lj*) (^--.ei)*^i|(^i>=i)' (4-11) 

(drnlj^k+i = {dTZl,^ri., + /.E^^(^\i)*(^^-.ei)*Gi^(Gi,Hi), (4.12) 

j=i ^ 



Ml 



di 



{Gl,El). 



(4.13) 



for i = 1, • • • ,s and k = 0, • • • , N — 1, if and only if it is a critical point of the 
function Sd : Cd{gi,g2) 1^; that is, dsd{cd) = 0. Moreover, the discrete flow map 
defined by the above scheme, : T^p Ihp, preserves the symplectic form wj^^ . 



Proof. Set r)k = g^ ^Sgk and Hf, = G^, SG\. The differential of Sd{cd) in the 
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direction z = {{bg^, 6fik}, {^Gl, J is given by: 

^''i-z= Y^Y.hhGl^{GlEl) . Hi + hb,-{GlEi) • SE^ 

k=0 i=l ^ ^ 



+ h(Sf^i,T-\g,'Gl)/h-Y,a 



'3^k 



l^k 



Collecting terms with the same variations and summation by parts using the bound- 
ary conditions 6 go = 6g]\f = gives, 

N-l si s \ 

k=l i=l \ j=l I 



hlSnk+uT ^{gk^gk+i)/h-'Yh 



h (^bi — {Gl, El) - ^ aj.fil - b^fik+i,dEl 

{drZl^.Tf^l + hb,Gl^^{GlEl),Hi 

+ (^-idr^ljyk+i + {dT:l^J^^k - EC^^^^^Q.OVi, 

Since dsd{cd) = if and only if dsd • 2 = for all z G Tc^Cd, one arrives at the 
desired equations with the elimination of and the introduction of the internal 
stage variables = dl / d^{G\^E],) for i = 1, • • • , s. Conversely, if q satisfies 
(4.9)-(4.13) then dsd{cd) = 0. 

Consider the subset of Cd given by solutions of (4.9)-(4.13). Let Sd denote the 
restriction of Sd to this space. Since each of these solutions is determined by a point 
in T/jp, one can identify this space with T/jp, and hence, Sd '■ Ihp ~^ Since Sd is 
restricted to solution space, 

dsd(ffo,^o,^o) • ibgo,5(,o,6no) = (^{dTZl^^)* hn, g^^bgN"^ - <;|(drr^^JVo,5'o ^'^fo) • 
Preservation of wj^^ follows from d^s^ = 0. ■ 
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The external and internal stages of (4.9)~(4.13) define update schemes on G x q* 
and (0 X X g*)*, respectively. 

Stormer-Verlet Integrators on Lie Groups. The generalization of the Stormer- 
Verlet method to Lie groups is given by evaluating (4.9)-(4.13) at the following 
two-stage RK tableau (implicit trapezoidal rule), 












1 


1/2 


1/2 




1/2 


1/2 



Given h and {gk,tJ-k), the method determines (3^+1,^^+1) by solving the following 
system of equations: 

M;/2 = |(9fc,Hi), (4.14) 
<' = f (9m,Hi), (4.15) 
(drnl.XMl/' = {drZl^Jf^k + ^9k^gi9k, ^1), (4.16) 
9k+i=9kT(h^{El + El)Y (4.17) 



^ik+i = mI'^ + -(dr_;,5,^J*5fc+i — (gfc+i,H2). (4.18) 
In particular, one uses the following procedure: 



• Eliminate gu+i in (4.15) using (4.17). Then solve for M^^^, S^, and using 
(4.14)-(4.16). This update is in general implicit. 

• Update gk+i using (4.17). This update is explicit. 

• Solve for ^k+i using (4.18). This update is explicit. 

Observe that if the Lagrangian is separable, then (4.16) is not implicit in the po- 
tential force term and one does not need to eliminate gk~\-i in (4-15) using (4.17). 

Variational Euler on Lie Groups. The variational Euler schemes come from 
evaluating (4.9)~(4.13) with the following tableaus: 









1 


1 




1 ' 


1 



forward Euler backward Euler 

The corresponding VPRK action sums take the following simple forms: 

= T.k=o h [iigk , Ck ) = T,k=o h [^{9k+i ,ik+i) 

+ {l^k,r~^{gk'^9k+i)/h - ^k)] + {fik+i,T~^{gk'^9k+i)/h - ^k+i) 
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Given h and (g^, the forward variational Euler method determines {gt+i, /"fe+i) 
by solving the following system of equations: 

9k+i = gkT{Hk), 
' i'^^hl+J^k+i = {dTZl^J*fik + hgk+i^{gk+i,^k+i), (4.19) 

fj-k+i = ^{gk+i,Ck+i)- 

The backward variational Euler method determines {gk+i, fJ-k+i) by solving the fol- 
lowing system of equations: 

9k+i = gkT{h^k+i), 
' (dTf^lj^k+i = {dTZl^J*f,k + hgklligk,Ck), (4.20) 
fJ-k+i = ||(5fc+i,&+i)- 

Euler-Poincare Integrators. In the case when the Lagrangian is G-left- invariant, 
the angular momentum updates in the above methods are identical and given by: 

(d^hLX^"'^^ = (drll^.y^^'^ (4-21) 

Examples. We now give various examples of Euler-Poincare integrators by making 
different choices of the map r and evaluating (4.21). 

(a) Matrix exponential. Suppose 

T = exp(0, t:q^G, 
which is a local diffeomorphism. 

Using standard convention the right trivialized tangent of the exponential map 
and its inverse are denoted by dexp : x g — > g and dexp^^ ^ ^ 0) ^-nd 
are explicitly given by, 

OO ^ OO 

dexp(x)y = ^ — — — ad^,y, dexp"^(x)y = ^ ^ ad^, y, (4.22) 
j=o ^ j=o ^' 

where Bj are the Bernoulli numbers] see §3.4 of Hairer, Lubich, and Wanner 
[2006] for a detailed exposition and derivation. 

Hence, (4.21) takes the form, 

{de^V'\Kk)T lik = (dexp-H-/iefc-i))Vfc-i- (4.23) 

(b) Pade (1,1) approximant. Suppose 

r{i) = cay(e) = (e - i/2)-\e + (4.24) 

which is the Pade (1,1) approximant to the matrix exponential and better 
known as the Cayley transform. The Cayley transform maps to the group for 
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quadratic Lie groups {SO{n), the symplectic group Sp{2n), the Lorentz group 
SO{3, 1)) and the special Euchdean group SE{3). 

The right-triviahzed tangent of the Cayley transform and its inverse are writ- 
ten below 

dcay(x)y = (e — x/2)~^y{e + x/2)~^, dcay~"^(x)y = (e — x/2)y{e + x/2). 

(4.25) 

For a derivation and exposition the reader is referred to §4.8.3 of Hairer, Lu- 
bich, and Wanner [2006]. Using these expressions (4.21) can be written as, 

h h 

+ (4.26) 

(c) Pade (1,0) or (0,1) approximant. Rather than use the exact matrix expo- 
nential one can use a Pade approximant, e.g., the Pade (1,0) approximant 

exp(^) « e -h C 

or Pade (0,1) approximant 

exp(0«(e-e)-^ 

However, since a Pade approximant is not guaranteed to lie on the group one 
needs to use a projector from GL{n) to G. In what follows G = SO{n) will be 
considered where a natural choice of projector is given by skew symmetrization. 

Suppose 



9-9^ 



2 

which comes from a first order approximant to the matrix exponential. This 
map is a local diffeomorphism from a neighborhood of e to a neighborhood 
of and its differential is the identity. Its right trivialized tangent can be 
computed from its derivative: 

By definition of the right trivialized tangent of r^-*^, it then follows that, 

dskew(x)(y) = ^^(^) - (4.27) 



Cardoso & Leite [2003] obtained the following theorem that explicitly deter- 
mines t(^). Moreover, they give necessary and sufficient conditions for its 
existence. 
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Theorem 4.10. Given £ so(n), a special orthogonal solution to the equation 

^ = 2 

can be written as 

1 /2 

where (^^ + e) is a symmetric square root. 

Proof. Since the skew-symmetric part of g is ^, one can write 5 as a sum of 
^ and a symmetric matrix S, 

r(e)=5 + 6 
Observe that ^ commutes with r(^) since 

2CTiO = (r(0 - miriO = r{Cf - e = 
Moreover, S satisfies an algebraic Riccati equation because, 

r{0*r{0 =e =^ + - - {^^ + e) = 0. 
And since ^ commutes with 5 (because it commutes with g), 

S' = {e + e), 

which completes the proof. ■ 
Hence, (4.21) can be written as, 

2 

2 

+ \ ad|^. ixk + ^ ad|^_^ (4.28) 

5 Conclusion 

In this paper a left-trivialized Hamilton-Pontryagin principle is derived for mechan- 
ical systems on a Lie group G. If the Lagrangian is left-invariant with respect 
to the action of G, it is shown that this left-trivialized HP principle unifies the 
Euler-Poincare and Lie-Poisson descriptions. In addition to its utility for implicit 
Lagrangian systems, the paper shows that this principle provides a practical way 
to design discrete Lagrangians. In particular, the paper explains how one can 
discretize the kinematic constraint using a Runge-Kutta Munthe-Kaas (RKMK) 
method. The paper shows that this leads to a novel generalization of variational 



5 Conclusion 22 



partitioned Runge-Kutta methods from flat spaces to Lie groups. In particular, one 
can generalize variational (or symplectic) Euler and Stormer-Verlet methods to Lie 
groups in this fashion. These methods inherit many of their attractive properties 
on flat spaces: efficiency, order of accuracy, symplecticity, symmetry, etc. 

Part II of this paper will develop a basic numerical analysis of these methods 
and report on numerical experiments on a class of nonreversible mechanical systems 
on Lie groups (moving rigid body systems) and chaotic dynamics of an underwater 
vehicle. To be specific the paper will: 

• prove order of accuracy of the VPRK integrators presented in this paper by 
invoking the variational proof of order of accuracy Marsden & West [2001]; 

• explain the numerics behind the Poincare sections provided in Figure 2.1; 

• demonstrate the superiority of these VPRK integrators compared to symmet- 
ric rigid body integrators when applied to a nonreversible system such as a 
rigid body on a turntable. 
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